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Abstract 

We study diffusion of particles in large-scale simulations of one-dimensional stochastic sandpiles, 
in both the restricted and unrestricted versions. The results indicate that the diffusion constant 
scales in the same manner as the activity density, so that it represents an alternative definition of 
an order parameter. The critical behavior of the unrestricted sandpile is very similar to that of its 
restricted counterpart, including the fact that a data collapse of the order parameter as a function 
of the particle density is only possible over a very narrow interval near the critical point. We also 
develop a series expansion, in inverse powers of the density, for the collective diffusion coefficient in 
a variant of the stochastic sandpile in which the toppling rate at a site with n particles is n(n — 1), 
and compare the theoretical prediction with simulation results. 
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I. INTRODUCTION 

Sandpile models are the prime example of self-organized criticality (SOC) , or scale- 
invariance in the apparent absence of control parameters In sandpiles, SOC arises 
via a control mechanism that forces the system, which possesses an absorbing-state phase 
transition, to its critical point 0, 0]. SOC in a slowly-driven sandpile corresponds to an 
absorbing-state phase transition in the conserved sandpile, which has the same local dynam- 
ics, but a fixed number of particles 

mm 

, [lO| . Conserved sandpiles are characterized 



by a nonconserved order parameter (the activity density 



field that does not evolve in regions devoid of activity llj. This class, known as conserved 



which is coupled to a conserved 



directed percolation (CDP), is distinct from that of standard directed percolation 12|. 

In recent years considerable progress has been made in characterizing the critical proper- 
ties of conserved stochastic sandpiles, although no complete, reliable theory is yet at hand. 
As is often the case in critical phenomena, theoretical understanding of scaling and univer- 
sality rests on the analysis of a continuum field theory or Langevin equation (a nonlinear 
stochastic partial differential equation) that reproduces the phase diagram and captures 
the fundamental symmetries and conservation laws of the system. Important steps in this 



direction are the recent numerical studies of a Langevin equation |12l . [iSj for CDP. The 
critical exponent values reported in Ref. jjj] are in good agreement with those found in 
simulations of conserved lattice gas (CLG) models 18|, ll9|, which exhibit the same symme- 
tries and conservation laws as stochastic sandpiles. The Langevin equation exponents are 
also consistent with the best available estimates for stochastic sandpiles in two dimensions 
[l^ . There is now good evidence that the one-dimensional stochastic sandpile belongs to 

the CDP universality class QE3. 

In this work we focus on an aspect of sandpiles that has received relatively little atten- 
tion: diffusion. Since the dynamics in these models involves hopping of particles between 
neighboring sites, one expects the particle diffusion constant D to follow a scaling behavior 
similar to that of the usual order parameter p. (A site is active if it bears two or more 
particles.) Here D is defined via the relation ((Ax)^) = 2Dt, where Ax is the particle dis- 
placement. We determine the scaling properties of the diffusion constant in extensive Monte 
Carlo simulations. 

Theoretical studies of the particle diffusion coefficient are hampered by the fact that po- 
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sitions of specific particles are not accessible in the usual stochastic description; instead the 
master equation describes the evolution of the probability distribution on the set of occupa- 
tion numbers {rij}. It is however possible to determine the collective diffusion coefficient 
by studying how a density perturbation Ap oc e*^^ relaxes. Dc is related to the relaxation 
time, Tfc, of this mode via r^, = 1/ {Drk^), in the small-/c limit. Using the path-integral based 
perturbation theory developed in 23|, we calculate the first three terms in the expansion of 
Dc in inverse powers of density p, for the stochastic sandpile in which the toppling rate is 
n{n — 1). 

The balance of this paper is organized as follows. In Sec. II we define the three models 
of interest. Sec III reports simulation results on D and p for two of these models. In Sec 
IV we develop a series expansion for the collective diffusion coefficient in a third model, and 
compare the predictions with simulation results. We close in Sec. V with a summary and 
discussion. 



II. MODELS 



We study three versions of the one-dimensional conserved stochastic sandpile, related to 



Manna's model 



20|. In these systems the configuration is defined by the set of occupation 



variables ni, ...,nL, giving the number of particles residing at each site on a ring of size L. 
All three versions are continuous-time Markov processes, in which transitions involve the 
"toppling" of an active site, i.e., one with rii > 2. The particular features distinguishing the 
three models are as follows. 



Basic unrestricted model (I) 21j. Each active site has a rate of unity to topple. When 
site i topples, two particles are transferred from this site to its neighboring sites {i — 1 and/or 
i + I, with periodic boundaries). The two particles jump independently; they jump to the 
left or to the right with equal probabilities. 



Restricted model (II) IJ, ll6|. The dynamics is that of model I except that no site may 



have more than two particles. If, when a site topples, a particle attempts to jump to a site 
bearing two particles, it returns to the toppling site. 

Modified unrestricted model (III). 23|. The dynamics is that of model I except that the 
rate at which a site topples is given by nj(?T,j — 1). 

While somewhat less convenient for simulation, model III is better suited to operator- 
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based theoretical approaches. Model II features a smaller set of states, rendering it more 
convenient for analysis via cluster approximations E^L-^ There is clear numerical evidence 
that model II belongs to the CDP universality class [14']; models I and III share the same 
symmetries and conserved quantities as model II and so are expected to belong to the CDP 
class as well. 

In conserved sandpiles, the particle density p serves as a temperaturelike control pa- 
rameter. Below a certain critical value, pc, the system eventually falls into an absorbing 
configuration (i.e., one devoid of active sites), while for p > p^, activity continues indefi- 
nitely, in the infinite-size limit. The order parameter associated with this absorbing-state 
phase transition is the activity density p, given by the stationary mean fraction of active 
sites in models I and II, and the stationary average {ni{ni — 1)) in model III. Numerical 
studies strongly support a continuous transition at pc, best estimates for the critical density 
Pc are 0.9488, 0.92978, and 0.9493, in models I, II, and II, respectively. (Note that for p > 1 
active sites always exist; pc, however, is strictly less than unity.) 




InA 



FIG. 1: (Color online) Asymptotic diffusion rate versus A = p — pc in model I, system sizes as 
indicated. Error bars are smaller than the symbols. 
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FIG. 2: (Color online) Asymptotic diffusion rate versus A in model II. 

III. SIMULATION RESULTS 

In this section we report simulation results for the particle diffusion rate D and the 
activity density p in models I and II. For each particle j, let Aj(t) = hjit) — hj{t), where 
h'j^lt) and hj{t) are the numbers of hops taken by particle in the positive (respectively, 
negative) directions up to time t. Then the particle diffusion rate is defined through the 
relation 



{[AM') = '^Dt (1) 

where the average is over all particles and (in principle) all histories of the system of = pL 
particles on L sites, starting from a given initial configuration or class of configurations. (In 
practice we generate initial configurations by adding particles randomly to the system, with 
the prohibition, in the case of model II, of triple or higher occupancy. We average over a 
set of Ng realizations of the process.) Note that the diffusion rate D defined above depends 
in general on the time t as well as on L and p. Determination of D in simulations requires 
that we store the particle displacements, which is not necessary if we merely wish to study 
the activity density. Since A{t) refers to a given history of the process, the quasistationary 
simulation method used in JL4| is not applicable here. 

Each time particle j hops, Aj changes by ±1, so that {[Ajit)]^) = {hj'(t) + hj(t)) = 
{hj{t)), i.e., the mean number of jumps up to time t. Particle j must be at an active site 
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in order to jump, but {hj{t)) is not simply equal to 2p, as would be the case if particle j 
were always to jump each time the site it resides at topples. In model I, for example, the 
probability of a given particle jumping is 2/n, if it is one of n particles at the toppling site. 
In model II, the particle always tries to hop when its host site topples, by it may be unable 
to move to the target site. Finally, in model III a particle at a site with occupation n hops 
at a rate of 2 (n — 1), so that the hopping rate actually grows with the occupation number n. 
Since the occupancies of nearby sites are correlated, the waiting times between successive 
displacements of a given particle are not independent. For these reasons, the relation between 
the hopping rate and the activity density involves subtle effects, different in each of the three 
models studied. It is nevertheless reasonable to expect that, as p — > Pc-, the scaling behavior 
of the diffusion rate will parallel that of the activity density. In particular, we might expect 
p and D to be governed by the same set of critical exponents at the transition. 

We simulate models I and II on rings of L = 6250, 12500, 25000, and 50000 sites, using 
eight independent realizations for the smallest size, six for L = 12500, and four for the two 
largest sizes. The studies are run for 10^ to 6 x 10^ time units, with the longest simulation 
times near the critical point. Each particle is assigned a label so that its cumulative dislo- 
cation Aj can be followed during the evolution. In model II, when two particles attempt to 
jump to the same site, and this site is singly occupied, one of the two particles is chosen at 
random to move to the target site, while the other remains where it is. 

We monitor the diffusion rate D{t), defined in Eq ([T]), and confirm that it approaches a 
stationary value at long times. Figs. [1] and [2] show the stationary value, D, as a function 
of A = p — pc, for models I and II, respectively. Several aspects of these results are worth 
commenting on. First, for the sizes considered here, D is apparently well converged to 
its limiting (L —>■ oo) value for A > 0.0025. Second, even for values of A such that the 
diffusion rate has converged, the slope of -D(A) on logarithmic scales changes appreciably 
with A, making a reliable estimate of f3 difficult. Finally, in model II, D drops sharply as p 
approaches 2: due to the height restriction, most particles cannot move. 

The corresponding results for the activity density p are shown in Figs. [3] and HJ respec- 
tively. The behaviors of D and p in both models appear quite similar, an impression that 
is confirmed in Fig. [5], which compares both quantities (in both models), for the largest 
system studied. Near the transition, p and D are virtually identical in the unrestricted 
model, while in model II they appear to be proportional. It is evident that neither p nor D 
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FIG. 3: (Color online) Stationary activity density p versus A in model I. 



can be characterized as following: a simple power law, an observation already made for the 



order parameter in model II in 



Note that for the system sizes studied here, there is no discernable finite-size effect for 
p — Pc ^ 0.0025. Since D{A) and p(A) do not follow simple power laws in this regime, 
there is no possibility of maintaining the data collapse under the usual kind of FSS scaling 
plot, that is, of p* = L^/'^-Lp versus A* = L^^'^-^A. As noted in IJ], a data collapse can 
only be achieved in the regime very near the critical point (i.e., A < 0.0025); and example 
of a data collapse for the diffusion rate data is shown in Fig. O In model I the data 
collapse is best using exponent values P = 0.289 and z/_l = 1.35; the corresponding values 
in model II are (3 = 0.285 and u± = 1.355. These values are consistent with those reported 
in [l^: f3 = 0.289(12) and z/^ = 1.355(18). We may therefore affirm, with a high degree of 
confidence, that models I and II belong to the same universality class, and that the diffusion 
rate and the order parameter exhibit the same critical scaling properties. 

In the paradigmatic examples of absorbing-state phase transitions, such as the contact 



process 25|, U^ , starting from a spatially homogeneous initial distribution, the order param- 
eter exhibits an initial power-law decay, p ~ t~^, at the critical point, before saturating at 
a quasistationary value. (Note that the power-law portion of the evolution is independent 
of system size.) It is of interest to know whether the order parameter and the diffusion rate 
exhibit similar behavior in the stochastic sandpile. Our results (Fig. [T]) for L = 50 000 show 
p and D decaying with an exponent 6 = 0.153(5). 
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FIG. 4: (Color online) Stationary activity versus A in model II. 



We also perform simulations of the spread of activity with time. In this case, a single site 
is given two particles initially, while the remaining N — 2 particles are distributed at random, 
one per site, over the rest of the lattice. In the contact process at criticality, starting with a 



single active site, the number of active sites grows as n{t) ~ |25l . 1261]. In the present case 
we find that both D and p follow an approximate power law with an exponent rj = 0.34(1), 
as shown in Fig. [HI The spreading exponents 6, rj, and Zgp are expected to satisfy the 
hyperscaling relation 45 + 2// = dzgp, which, using Zsp = 2/z, with z the usual dynamic 
exponent, can be written as rj = 1/z — 26. Using the value of 6 cited above, and z = 1.50(4) 
from Ref. [l^, this yields rj = 0.36(3), which is consistent with our numerical estimate. 

IV. COLLECTIVE DIFFUSION COEFFICIENT: THEORY AND SIMULATION 



In this section we apply the operator formalism and perturbation theory derived in [23 1 
to the evaluation of the collective diffusion coefficient of model III on a ring of N sites. 
We begin by writing the master equation for the process in the form 

where 

l^) = Ep(M'i)IW) (3) 

{n} 
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FIG. 5: (Color online) Stationary activity (open symbols) and asymptotic diffusion rate (filled 
symbols) versus A in models I (squares) and II (circles), system size L = 50 000. 

is the probability distribution. Here p({n},t) is the probability of configuration {ra}, and 
the state |{n}) is a direct product of states |nj), representing exactly rij particles at site j. 
These states are normalized so: {n'\n) = n\6n,n'- 

Defining creation and annihilation operators via the relations, 

ailm) = ni\ni-l) (4) 

and 

TTi\ni) = \ni + l), (5) 
so that [aj,7rj] = 6ij, the evolution operator for the one-dimensional stochastic sandpile is 

^ = E[^(^-i + ^mr-^'l4 (6) 

i 

Since the system is translation-invariant it is convenient to introduce the discrete Fourier 
transform via 
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In A* 



FIG. 6: (Color online) Scaled diffusion rate D* = Lf^/^^D versus scaled distance from critical point 
A* = Li/'^iA in model I. 




FIG. 7: (Color online) Initial decay of activity (upper) and diffusion rate (lower) at criticality in 
model I; system size L = 50 000. The slope of the straight line is -0.153 



«fc = E e '^^Cij , (7) 
j 

with inverse 

«. = ^Ee^^V, (8) 

k 
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FIG. 8: (Color online) Initial growth of activity (squares) and diffusion rate (circles) at criticality 
in model II, starting with a single active site; system size L = 50 000. 

(and similarly for other variables), where the allowed values of the wavevector are: 

27r 27r 2tx 2tx . . 

A. = -vr, -. + -,...--, 0, (9) 

(To avoid heavy notation, we indicate the Fourier transform by the subscript k; the subscript 
j denotes the corresponding variable on the lattice.) In the Fourier representation, the 
evolution operator takes the form 



L = -N ^ ^ UJki^k2T^kiT^k2(^k3a-ki-k2-k3, (10) 
ki,k2,k:j 



where ujki,k2 = 1 ~ cos /ci cos /i;2. As explained in Ref. [23|, the evolution operator may be 
rewritten as 



with 



and 



L = Lo + L, (11) 



Lo = ~N~^Y.^k7r-kak, (12) 

ky^O 



^k-i,,—k\—k2—kz ^fea fci — ^2 — ^3 ^k\ ^k2 

ki,k2,k3¥=0 

- 2pN~'^ ^ t^fc2,-fci-fc2^fc2^-fcl-fc2'^fcl 

ki,k2¥=0 

- 2N~^ LLi_ki-k2,o'^-ki-k20'kiak2 - P^N~^J2^k-k'n:k'n:-k, (13) 

fei,fc27^0 k^O 
11 



where 



7fc = 4:pujkfl = 4p(l - cos k). (14) 

This transformation is based on the observation that, due to particle conservation, the 
operator A^~^ Y.j be equated to the particle density p. In Eq. f lTSl) . it is understood 

that none of the wavevectors associated with the operators a and vr may be zero. 

Let Pn = e~^p"/n! denote the Poisson distribution with intensity p, and define \P)i = 
J2nPn\n)i as the Poisson-distributed state at site i. Then the uniform product-Poisson 
distribution is \P) = (S)j|P)j. The latter is an eigenstate of the diffusion operator with 
eigenvalue zero, i.e., "Pj-P) = 0, where 

P = ^ ^ [vTj-i - 2TTj + TCj+i] aj = -A^~^ ^kflT^-kak, (15) 
represents nearest-neighbor hopping at unit rate. (Note that Lq = 4pP.) 



To study collective diffusion, we consider an initial condition in which the uniform 
Poisson-product is weakly perturbed by a density modulation with wavevector k: 



1^(0)) = N-''Kk\P) (16) 



Introducing the notation. 



(i-En^^Ki (17) 

{n} j "-i- 

for the projection onto all possible configurations, the mean number of particles at site j is 
given by 

Ut) = {nM = ( k.l^W) (18) 
or equivalently, in the Fourier representation, 

Mt) = ( |afe|^(t)). (19) 
Note that for the initial distribution of Eq. [161 with 7^ 0, we have 
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FIG. 9: (Color online) Projection <I>(t) in simulations of model III. System sizes (lower to upper) 
N = 400, 800,..., 6400, particle density p = 3.0. The slope of the straight line is -1/2. 

N-'{ IsvTfclP) = (20) 

where we used the relations [ag,7rfc] = N6g^^k and ( \iTk\P) = N6kfi- Thus (p-kif) represents 
the amphtude, at time t, of the density perturbation created at time zero. 

We assume that for long times and long wavelengths the mean density (f)j satisfies the 
diffusion equation d(f)j/dt = DcA'^(f)j with the discrete Laplacian, leading, in the small-fc 
limit, to (pkit) — 0fc(O) exp[— Defeat]. Letting (pki^) denote the Laplace transform, we have, 
in the small- 2; limit, (j)k{z) — l/{z + Dck"^), so that, 

= lim ^r^— -. (21) 
Laplace transforming the formal solution of the master equation we find 
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FIG. 10: (Color online) Projection <I>(t) in simulations of model III with = 12800 and p = 
The solid line is given by Eq. (j37p with Dc = 3.76. 



(P_k = N-^{ \a.k^-Y^k\P) 
z — L 

We may develop the solution in a series in powers of 1/p using the operator identity 



1 1 1^1 1 ^ 1 ^ 1 

+ — + — Li — Li — + 



z — Lq — L\ Z — Lq z — Lq Z — Lq z — Lq z — Lq Z — Lq 

Evaluation of the contributions to this series is facilitated by use of the identities Lq\P) 
ak\P) = NpSk,o\P), and, 

( laqTTfclP) = N6g-k + N'^pSqfihfi. 

Note also that Lott^Ip) = —'~fk'^k\P), and in general, 

LoUk^ ■ ■■TikjP) = -Suk^ ■ ■ -vrfc^lP). 
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where 5 = + ■ ■ ■ + 7a;„ . Thus Lq may be inverted on the space of states of the form 
TTfei ■ • • vrfc,jP) provided that not all of the wave vectors are zero: on this space Lq ^ is simply 
—1/5 times the identity operator. 

The first term in the expansion of Eq. fl2^ is readily evaluated as, 

= N-\ |a., n,\P) = — ^, (26) 



which gives limfc^^^o k'^4>-k{z) = l/{2p) + Oil/p^). Subsequent terms in the expansion 



may be evaluated using the diagrammatic perturbation approach developed in [23|]. In this 
representation each term in Li corresponds to a vertex, with operators corresponding to 
lines entering the vertex at the right, and operators vr^ to lines leaving at the left. Each line 
that leaves a vertex must be joined ("contracted") with a line entering some other vertex to 
the left. The operator Li, Eq. (fT3|) . consists of four parts or vertices, designated respectively 
as a crossing (two lines in, two out), a bifurcation (one in, two out), a conjunction, and a 
source. We denote these contributions as L^, L;,, Lc and Ld, respectively. There are two 
diagrams that contribute to (f)-k{z) at order 1/p^. One arises from the term. 



,r1/l Ir Ir ^ Ir^K 1 — COS ^ COS ( — g) 



z — Lq z — Lq z — Lq 8A''p2 ^ 2 — cosg — cos(A; — g) 



1 



16p^(l — COS /c) 



(27) 



while the second is, 



N-\ \a_,^-L,^-L,^-Ld^-n,\P) = -—— -. (28) 

z — Lq z — Lq z — Lq z — Lq 32p^(l — cos /ej 



At order there are 17 diagrams, leading to 



-^c 1 _i_ 1 I 0.1088899 I (^^) 
+ 8p + p2 
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In 23| the stationary activity density p = — 1)) was found to grow asymptotically 
as p^, with correction in inverse powers of p; here we find that Dc grows only linearly with 
p. For comparison we write the results for p and D^. in the form: 



P = P 



1 0.0492525 
1-- ^ + 



4j9 



and 



1 0.093265 
1-:^ + 



5p 



(30) 



(31) 



These expressions are reliable for 4p ^ 1 but cannot of course be applied in the vicinity of 
the critical density, pc — 0.9493. 




FIG. 11: (Color online) Collective diffusion constant versus particle density p in model III. 
Squares: simulation; line: theory, Eq. (j29p . The crosses denote simulation values for the stationary- 
state collective diffusion constant Z)c,s- 
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A. Comparison with simulation 



We determine the collective diffusion coefficient in model III via analysis of the projection 
of the configuration {n{j,t)} on the initial configuration. Recalling that {nj{t)) = p, the 
particle density, we let fj(t) = rijit) —p denote the excess particle number at site j and time 
t. Initially, the {uj} are independent, Poisson-distributed with mean p. Consider now 

(my) ^ ^ 

where the angular brackets denote an average over sites and over realizations, including the 
random initial configuration. For the set of wavevectors k defined in Eq. ([9]), let ipk(t) 
denote the discrete Fourier transform, 

N 

Mt) = j:fAty"' (33) 



In the small-A; limit we expect ipk to follow, 

Mt) = MO)e-''^'"' (34) 

Using the fact that the fj{0) are independent, zero-mean random variables with var[/j(0)] = 
p, it is straightforward to show that 

m = 4 E e"""^"'* r e-^=''*rfA; (35) 

N ^ ZTT J-TT 

For times such that Dj^'^t ^ 1 we may extend the limits of integration to ±oo, yielding 



Thus if $(t), as determined via simulation, can be fit for large t with an expression of the 
form A/t^^"^, then Dc = l/^AnA^). In practice, however, a more reliable procedure is to fit 
the simulation data to the full lattice expression 



17 



which involves the single adjustable parameter Dc, and is capable of fitting the data at short 
as well as long times, and for various system sizes. 

We determine on rings of = 200,400,800, ...,25600 sites, for particle densities p 
in the range of 1 to 3. Fig. [H] shows that as the system size increases, approaches a 
power-law decay with an exponent of 1/2; an example of data fit by Eq. (1371) is shown in 
Fig. [TOl The estimates for Dc, obtained by fitting J-'{t) to the simulation data, are compared 
with the theoretical prediction, Eq. fl29|) . in Fig.[TTl the agreement is quite good for densities 
p > 2. In Fig. [m the relatively large error bars associated with the simulation results for 
p = 1.1 and 1.2 reflect the fact that the simulation results for are less well fit by the 
theoretical expression, Eq. (1571) . than for other particle densities. Curiously, for p = 1, 
despite being nearer the critical point, the fit is again quite good. 

A similar analysis can be applied to extract the value of Dc,s from simulations in the 
stationary state. In this case, we allow the system to relax, so that the configuration at 
time zero is typical of the stationary distribution. Now, however, the /j(0) are no longer 
independent, Poisson distributed variables, and the power spectrum of fluctuations (1^3^(0)^) 
is no longer constant. We therefore fit the data for using the expression 



with (|0fc(O)p) determined via simulation. The resulting stationary values of are close 
to, but slightly greater than, those found using the Poisson initial distribution (see Fig. [TTl) . 
It is worth noting that in the stationary state the projection $(t) appears to decay with 
a power smaller than 1/2 (a typical exponent value is about 0.41). This does not imply 
anomalous behavior as the data can again be fit using the hypothesis fk(t) = ipk{0)e~^''^^^. 
For densities p < 1.2 however, the simulation data are not well fit by the function J-'s{t). In 
this regime the Fourier amplitudes {(p-k{t)(pk{0)) (calculated in simulations) do not follow 
a simple exponential decay. (The data suggest a crossover to stretched-exponential decay 
at long times.) Thus, near the critical point, we find evidence of anomalous relaxation, as 
previously noted in stochastic sandpiles 
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V. SUMMARY 



We study diffusion in stochastic sandpiles. In the first part of this work we determine 
the particle diffusion coefficient in sandpiles in which all active sites share the same toppling 
rate. We find, in both the restricted and unrestricted cases, that the diffusion constant scales 
in the same manner as the order parameter (the activity density). Our results confirm that 
the restricted and unrestricted models belong to the same universality class, and that both 
models exhibit a finite-size scaling collapse of data over an unusually narrow region of the 
control parameter (that is, the particle density p). 

The second part of this study deals with a sandpile in which the toppling rate at site i 
is ni{ni — 1). In this case it is possible to derive a short series for the collective diffusion 
constant, starting from a Poisson-product initial state. The resulting expression compares 
well with simulation for densities well above Pc- The collective diffusion constant Dc is 
extracted from simulations using the projection of density fluctuations at time t onto their 
initial values. We expect this approach to be useful in determining Dc in other systems, 
such as interacting lattice gases. We defer a detailed investigation of collective diffusion in 
the critical region to future work. 
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